Atomic Layering at the Liquid Silicon Surface: a First-Principles Simulation 
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We simulate the liquid silicon surface with first-principles 
molecular dynamics in a slab geometry. We find that the 
atom-density profile presents a pronounced layering, similar 
to those observed in low-temperature liquid metals like Ga 
and Hg. The depth-dependent pair correlation function shows 
that the effect originates from directional bonding of Si atoms 
at the surface, and propagates into the bulk. The layering has 
no major effects in the electronic and dynamical properties of 
the system, that are very similar to those of bulk liquid Si. 
To our knowledge, this is the first study of a liquid surface by 
first-principles molecular dynamics. 

PACS numbers: 68.f0.-m, 61.25.Mv, 61.20.Ja, 71.22.+i 



Liquid metal surfaces have attracted much attention 
during the last years p]-f|. Their particular properties, 
very different from those of non-metals, have been inves- 
tigated by Angstrom-resolution experiments, and simu- 
lated by different approaches. One of their most inter- 
esting features is the atomic layering, a density oscilla- 
tion that originates at the sharp liquid-vapor interface, 
and extends several atomic diameters into the bulk. Al- 
though some experiments supported an increased surface 
density at /-Hg Q , and some kind of atomic layering was 
predicted theoretically , its existence has been demon- 
strated unambiguously only recently by X-Ray reflectiv- 
ity in Hg §] , Ga |,| and Ga-In alloys §. Most of 
the experiments have been done close to the low melting 
temperatures of these metals, but Regan et al Q studied 
the Ga surface up to 170°C. They found that capilarity 
waves strongly decrease the reflectivity peak heights, but 
not the peak widths, suggesting that the decay length 
of local layering is temperature-independent, and that 
surface layering is still present quite above the melting 
point. 

Different explanations have been proposed to account 
for this effect 
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Rice et al |5j] have argued that the 
abrupt decay of the delocalized electron density forms 
a flat potential barrier against which the ions lay or- 
derly, like hard spheres against a hard wall. Tosatti et 
al H have used the glue model of metallic cohesion to 
argue that surface atoms, trying to effectively recover 
their optimal coordination, alternatively increase and de- 
crease their density. Surface layering effects, like surface- 
enhanced smectic ordering, have also been observed in 



liquid crystals pT| . In this case, its origin is the ten- 
dency of the highly nonspherical molecules to present a 
particular orientation towards the surface. 

Liquid silicon (Z-Si) is a rather peculiar system. Sili- 
con transforms, at 1684 K, from a covalent semiconduc- 
tor solid, with diamond structure and coordination 4, 
to a liquid metal. Experiments JlJ,[l3| and MD simula- 
tions @,^5) show that its coordination (^6 — 7) is lower 
than that of typical liquids (~ 12), due to the persis- 
tence of directional bonding in the liquid phase. In spite 
of the enormous literature on its solid surfaces, very lit- 
tle is known on the structure of the liquid silicon surface. 
Measurements are very difficult because of its high reac- 
tivity and melting temperature. Model calculations, and 
computer simulations with semiempirical potentials, are 
also difficult because of the mentioned coexistence of co- 
valent and metallic bonding, and its unknown interplay 
at the surface. 

In this letter we present a study of the Z-Si surface 
by first-principles molecular dynamics (MD) simulation 
p6| . This approach deals equally well with covalent and 
metallic bonding, and it is therefore very well suited 
for this problem. Electrons are treated by solving the 
Kohn-Sham |l7]] equations selfconsistently for each ionic 
configuration, using the local density approximation for 
exchange and correlation. The quantum mechanically 
obtained forces are then used to generate the classical 
trajectories of the ion cores. The calculations were per- 
formed with the SIESTA program using a linear com- 
bination of numerical atomic orbitals as the basis set, and 
norm-conserving pseudopotentials p9|] . A uniform mesh 
with a planewave cutoff of 40 Ry is used to represent 
the electron density, the local part of the pseudopoten- 
tial, and the Hartree and exchange-correlation potentials. 
Only the T fc-point was used in the simulations, since pre- 
vious work [^0| found cell-size effects to be small. For the 
present calculation we used a minimal basis set of four 
orbitals (1 s and 3 p) for each Si atom, with a cutoff 
radius of 2.65 A. We have extensively checked the basis 
with static calculations of different crystalline Si phases 
and solid surfaces, and MD simulations of the bulk liquid 
pT| . The energy differences between solid phases are de- 
scribed within 0.1 eV of other ab-initio calculations. The 
diamond structure has the lowest energy, with a lattice 
parameter of 5.46 A (0.5 % larger than the experimental 
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value). Adatom- and dimer-based (111) and (100) sur- 
face reconstructions found in other ab-initio calculations 
p2[ are well reproduced, with geometries and relative en- 
ergies changing less than ~ 0.1 A and ^0.15 eV when 
moving from a T-point calculation with a minimal basis 
set, to a converged fc-sampling with double-C and polar- 
ization orbitals. The structural, electronic and dynamical 
properties of /-Si are in good agreement with other ab- 
initio calculations at the same density and temperature 
Hp). The calculated diffusion constant (1.5 x 10~ 4 
cm~/s) is somewhat smaller than that obtained with a 
double-^ or polarized basis (1.7 — 2. Ox 10 -4 cm 2 /s) which 
is in agreement with other ab-initio simulations. We in- 
terpret that the minimal basis overestimates the energies 
of the saddle point configurations occurring during dif- 
fusion, but we consider that this is not critical for the 
present application. Also, we leave for a future work the 
inclusion of spin fluctuations, which affect significantly 
the diffusion constant but not the structural properties 
of the liquid f|J. 

We first perform a long simulation of bulk /-Si at 
T=1800 K ^3|, using a cubic 64-atom cell with periodic 
boundary conditions. The fixed cell size (10.58 A) was 
adjusted to obtain zero mean pressure, and corresponds 
to a density 3% smaller than the experimental density 
near the melting point. We then construct our initial 96- 
atom slab by repeating one bulk unit cell in the x and y 
directions, and one and a half cells in the z direction, plus 
10 A of vacuum. No particles leave the slab during the 
30 ps simulation. After a relaxation of 10 ps, the system 
reaches equilibrium and the averaged magnitudes are es- 
sentially the same for the next or the last 10 ps, and for 
both sides of the slab. These long relaxation and obser- 
vation times are required because the calculated density 
autocorrelation time at the surface (~ 1 ps) is consider- 
ably longer than the typical bulk-liquid correlation times 
(~ 0.1 ps) 0. The average surface energy (836 ± 40 
dyn/cm) is in good agreement with the experimental sur- 
face tension (850 dyn/cm at 1800 K) |m]| , suggesting a 
small cntropic contribution. 

In Fig. @ (solid line) we present the ionic density pro- 
file p(z). It shows a pronounced atomic layering, with 
similar features as those reported for the Ga surface j|. 
Like in that case, p(z) can be fitted accurately by a sharp 
error function at the surface, and a sinusoidal wave with 
an exponential decay towards the bulk: by superimpos- 
ing two such functions (not shown), for both surfaces, we 
obtain similar values of the parameters and an oscillation 
period of 2.5 A. To check that the observed layering is not 
a sign of incipient crystallization, we have computed sev- 
eral magnitudes in the central region of the slab (\z\ < 3.7 
A). The radial and angular distribution functions, elec- 
tronic and vibrational densities of states, and the diffu- 
sion constant, are all very similar to those of bulk Z-Si, 
and have no resemblance of those in the solid phases. As 
an example, we compare in Fig. the bond-angle distri- 



bution function for the bulk and slab simulations, using 
a bond cutoff distance of r m =3.10 A p5j. For a better 
understanding of the origin of the layering, we compute 
the normalized density-density correlation function [E6| : 



c p {z ,z) = 



(Sp(z ,t)Sp(z,t)) 
(5pf*(zo,t))l(5p*(z,t))i 



(1) 



where () denotes time average and Sp is the difference 
between the instantaneous density at time t, p(r, t) = 
J2iLi <5(r — i"i (£)), and the average density p(r) — (p(r, t)), 
where 5(r) is Dirac's function. Fig. |l| also shows c p (zoj z ) 
for zq at the positions of the outermost peaks of each 
side. Its decaying oscillation is clearer than that of the 
density profile, all whose relevant features match very 
well with the superposition of the two c p 's. The appar- 
ent lack of decay of p(z) towards the interior is peculiar 
to our particular slab thickness, because the superposi- 
tion is positive at the center of the slab, and negative at 
the two surfaces. Most important is, however, that the 
two surface-induced oscillations are clearly independent 
of each other (c p 's out of phase), and incommensurate to 
the slab thickness. The density layering is thus an intrin- 
sic surface feature and not a result of finite size effects. 

In order to obtain information about the bond orienta- 
tions at the surface we calculate the two-particle density: 
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To represent /?2( r o> r )> we first average over the directions 
parallel to the surface: 



p 2 (z ;z,x)= 
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d 3 r' d 3 r'p 2 (r' ;r')6(z -z' ) 



5(z - z')S(x - - x' o r + (y> - y' ) 2 ), (3) 

where A is the area of the simulation cell. In Fig. || 
we show P2(zq; z, x) for zq located at the three peaks of 
p(z). Fig. ||(a) shows a clear tendency of surface atoms 
to form bonds parallel and normal to the surface. The 
height of the correlation peaks goes well beyond those of 
the density, which can be seen also in the figure (notice 
that p2(ro; r ) ~ * p( r ) for l r ~ r o| — * oo). This shows 
that the bond-induced correlations are responsible for the 
layering of the density, and not the other way around. 
Fig. |^(b) shows a similar, but attenuated tendency, that 
disappears in the third layer (Fig. ||(c)), which already 
has a very symmetric, bulk-like pair correlation function. 

Further insight can be obtained from the z-dependent 
coordination n(z), defined as the average number of 
neighbors within a distance r m . At the bulk, we ob- 
tain n=6.4, in agreement with the experimental value 
pjl . The distribution of local coordinations (DLC) is also 
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very close to those of other ab-initio calculations |1J,|15|, 
showing a maximum at coordination 6. We can also use 
the bulk simulation to construct an ideally terminated 
surface, cutting abruptly the system at, say, z—0. We 
then find n(z) = 4.3 at z — 1.0 A, which is the distance 
between the outermost peak and the inflection point in 
the slab density profile. At the outermost peaks of the 
actual slab, we obtain n(z)=5.3, and a DLC peaked at 
5. These values show that surface structural rearrange- 
ments increase the coordination of the ideally terminated 
surface, reaching a value of only one neighbor less than 
in the bulk. If we associate coordination 6, in the bulk 
liquid, with an octahedral arrangement, a simple picture 
can be drawn, in which the surface atoms try to preserve 
their bulk environment while minimizing the number of 
broken bonds. As a consequence, the octahedra get ori- 
ented in the surface so that only one broken 'bond' points 
towards vacuum, with another bond towards the interior 
and four bonds laying on the surface plane. This picture 
is consistent with figures ||(a) and ||(d). In the latter, we 
have restricted the sum over i in eq. (||) to particles hav- 
ing coordination 5, what results in even more pronounced 
peaks in the x and z directions. Also, we note that the 
maximum of /^(^o; 0) occurs at \z — zo\ = 2.5 A, what 
explains the same period observed in the density profile. 
The same distance is found for the in-plane surface bonds 
(i.e. for P2(zq; Zq,x)), and for the bulk bonds. Thus, con- 
trary to other metals, we do not find a shortening of 
the surface bonds, and the silicon surface layering seems 
to be related only to the bond orientations. However, 
it must be emphasized that the bond angle distribution 
(Fig. |^) is very wide, indicating a large variety of fluc- 
tuating atomic environments jl4j], so that our 'oriented 
octahedra' should be considered only as a very rough and 
qualitative picture. 

An interesting question is whether the surface struc- 
tural rearrangements produce a noticeable signature in 
the electronic structure. In Fig. |i[ we compare the lo- 
cal density of states (LDOS) at the outermost peaks of 
p(z) and at the center of the slab. Although fc-sampling 
is important for a converged LDOS |2l],|27j, we use here 
only the T-point eigenvalues, to facilitate the compari- 
son with previous work fl4|]l5|| and because we focus on 
its spatial variation. It can be seen that, apart from a 
slight narrowing, due to the diminished surface coordina- 
tion, there are no major differences, suggesting that the 
surface and bulk atomic environments are rather similar, 
and again pointing towards bond orientations as respon- 
sible for surface layering. 

In conclusion, we have performed the first study of a 
liquid surface by first-principles MD simulation. In spite 
of the high melting temperature of Si, we find a marked 
layering of the density near the surface, similar to those 
observed in other metals, like Ga and Hg, with low melt- 
ing temperatures. However, the surface layering of Si 
seems to have an origin at least partially different from 



that in other metals, with remanent directional covalent 
bonding playing an essential role. In spite of the rather 
slow decay of the layering towards the bulk, the average 
structural, dynamical, and electronic properties converge 
very rapidly to their bulk liquid values. Although more 
converged simulations would be highly desirable in the 
future, we consider that this work provides a new qual- 
itative understanding of the complex structure of liquid 
surfaces. 

We acknowledge useful discussions with E. Chacon and 
M. Weissmann. This work was supported by Argentina's 
CONICET and by Spain's DGES grant PB-0202. 
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FIG. 1. Solid line: density profile p(z) (relative to the bulk 
density po) of a liquid-Si slab, averaged during 20 ps and 
smeared by a gaussian convolution of 0.15 A. Dashed (z =6.2 
A) and dotted (zo=-6.2 A) lines: density-density correlation 
function c p (zo,z) (eq. (1)). We use a centered window of 0.5 
A for z , and a gaussian smearing of 0.5 A for z — z . 

FIG. 2. Distribution of bond-angles in the central region of 
the slab (continuos line) and in bulk /-Si (dashed line). 

FIG. 3. Two-particle density p2(zo; z,x) (eq. (3)) for: (a) 
2 = 6.2 A; (b) z = 4.0 A; and (c) z = 1.3 A. The 'volcano' 
is centered at x = 0, z — zo (position of the reference particle). 
p2 has been extended symmetrically to x < to facilitate its 
visualization (a line of ripples is produced by noise due to 
poorer statistics in eq. (3) at x ~ 0). (d) The same as (a), 
but restricted to atoms at zo having coordination 5. 

FIG. 4. Local density of electron states of the atoms in 
the outermost density peaks (solid line), and of those at the 
center of the slab (dashed line). 
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